function [N_D,delta_N_D,n_d]=newton_direction(u,delta_u)
global c a_t b p e t n_t n_p x_t y_t ar nth0 pwr;

[fu,dfu,sfu]=nonlinear_terms(u);

c1=1/c;
a=a_t-dfu;
a=a';
rhs=-c1*delta_u+a_t.*u-fu;
rhs=rhs';
n_d=assempde(b,p,e,t,c1,a,rhs);  %Solve -c\delta n_d+(a_t-dfu)n_d=rhs 
N_D=pdeintrp(p,t,n_d);          %for Newton direction at p and t.
delta_N_D=c*(-rhs+a.*N_D);
